rm(list=ls())

setwd("C:/Users/torewig_adm/Dropbox/Samarbeid/CHK, TW og SD/Opposition movements/FINAL SUBMISSION JOP/Replication materials")

library(foreign)
library(readstata13)
data <- read.dta13("ForCDE.dta")

names(data)

#subsetting to get complete cases
data <- subset(data, select=c(v2x_polyarchy_t1, urban_dominate2, indwork_dominate2, regviol, v2x_polyarchy, e_migdppcln,
                              e_mipopulaln, e_miurbani, nonviol, defect, regaid, viol))

data <- data[complete.cases(data), ]

#####CONDITIONING ON REGIME VIOLENCE

first <- lm(v2x_polyarchy_t1 ~ urban_dominate2 + indwork_dominate2 + regviol + v2x_polyarchy + e_migdppcln + 
              
              e_mipopulaln +  e_miurbani ,
            data = data)
summary(first)

data$direct <- data$v2x_polyarchy_t1-coef(first)["regviol"]*data$regviol

direct <- lm(direct ~ urban_dominate2 + indwork_dominate2 + v2x_polyarchy  + e_migdppcln + 
               
               e_mipopulaln +  e_miurbani , data = data)
summary(direct)


#bootstrap the standard errors
boots <- 1000
fl.boots <- rep(NA, times = boots)
for (b in 1:boots) {
  data.star <- data[sample(1:nrow(data), replace = TRUE),]
  boot.first <- lm(v2x_polyarchy_t1 ~ urban_dominate2 + indwork_dominate2 + regviol + v2x_polyarchy + e_migdppcln + 
                     
                     e_mipopulaln +  e_miurbani
    
    
    
    , data = data.star)
  boot.direct <- lm(direct ~ urban_dominate2 + indwork_dominate2 + v2x_polyarchy + e_migdppcln + 
                      
                      e_mipopulaln +  e_miurbani , data = data.star)
  fl.boots[b] <- coef(boot.direct)["urban_dominate2"]
}
sd(fl.boots)


boots <- 1000
fl.boots2 <- rep(NA, times = boots)
for (b in 1:boots) {
  data.star <- data[sample(1:nrow(data), replace = TRUE),]
  boot.first <- lm(v2x_polyarchy_t1 ~ urban_dominate2 + indwork_dominate2 + regviol + v2x_polyarchy + e_migdppcln + 
                     
                     e_mipopulaln +  e_miurbani
                   
                   
                   
                   , data = data.star)
  boot.direct <- lm(direct ~ urban_dominate2 + indwork_dominate2 + v2x_polyarchy + e_migdppcln + 
                      
                      e_mipopulaln +  e_miurbani , data = data.star)
  fl.boots2[b] <- coef(boot.direct)["indwork_dominate2"]
}
sd(fl.boots)


###########################################
#######       COEFFICIENT PLOT     ########
###########################################

est <- rbind(coef(direct)["urban_dominate2"],coef(direct)["indwork_dominate2"], coef(first)["urban_dominate2"], coef(first)["indwork_dominate2"] )

#conventional standard errors
ses <- summary(first)[["coefficients"]][,2]

se <- rbind(sd(fl.boots), sd(fl.boots2), ses[2], ses[3])


#plot
pdf("CDEfig.pdf", height=6, width=7)
plot(1, type="n", xlab="", ylab="", xlim=c(-.08, .08), ylim=c(0,10), axes=F, main="Effects conditional on regime violence")

legend("topleft", pch=c(19), cex=c(1), legend=c("CDE: IW dominate", "CDE: UM dominate", "Conventional coefficient: UM dominate", "Conventional coeffient: IW dominate"), 
       col=c("red", "green", "blue", "purple")
       
       )
abline(v = 0, lty = 3, col = "black") # add verticle line

#CDE's
lines(c(est[1,] + 1*se[1,], est[1,] - 1*se[1,]), c(4,4), cex=1.2, col="orange")  
lines(c(est[2,] + 1*se[2,], est[2,] - 1*se[2,]), c(5,5), cex=1.2, col="orange")  
lines(c(est[3,] + 1*se[3,], est[3,] - 1*se[3,]), c(6,6), cex=1.2, col="orange")  
lines(c(est[4,] + 1*se[4,], est[4,] - 1*se[4,]), c(7,7), cex=1.2, col="orange")  

points(est, c(4,5,6,7), col=c("red", "green", "blue", "purple"), cex=1.3, pch=19)

# add axes and labels                                                                                         # add bottom axis
axis(side= 1)
mtext(side = 1, "Coefficient", line = 3)      

dev.off()




#####CONDITIONING ON foreign support

first <- lm(v2x_polyarchy_t1 ~ urban_dominate2 + indwork_dominate2 + regaid + v2x_polyarchy + e_migdppcln + 
              
              e_mipopulaln +  e_miurbani ,
            data = data)
summary(first)

data$direct <- data$v2x_polyarchy_t1-coef(first)["regaid"]*data$regaid
summary(data$direct2)



direct <- lm(direct ~ urban_dominate2 + indwork_dominate2 + v2x_polyarchy  + e_migdppcln + 
               
               e_mipopulaln +  e_miurbani , data = data)
summary(direct)


#bootstrap the standard errors
boots <- 1000
fl.boots <- rep(NA, times = boots)
for (b in 1:boots) {
  data.star <- data[sample(1:nrow(data), replace = TRUE),]
  boot.first <- lm(v2x_polyarchy_t1 ~ urban_dominate2 + indwork_dominate2 + regaid + v2x_polyarchy + e_migdppcln + 
                     
                     e_mipopulaln +  e_miurbani
                   
                   
                   
                   , data = data.star)
  boot.direct <- lm(direct ~ urban_dominate2 + indwork_dominate2 + v2x_polyarchy + e_migdppcln + 
                      
                      e_mipopulaln +  e_miurbani , data = data.star)
  fl.boots[b] <- coef(boot.direct)["urban_dominate2"]
}
sd(fl.boots)


boots <- 1000
fl.boots2 <- rep(NA, times = boots)
for (b in 1:boots) {
  data.star <- data[sample(1:nrow(data), replace = TRUE),]
  boot.first <- lm(v2x_polyarchy_t1 ~ urban_dominate2 + indwork_dominate2 + regaid + v2x_polyarchy + e_migdppcln + 
                     
                     e_mipopulaln +  e_miurbani
                   
                   
                   
                   , data = data.star)
  boot.direct <- lm(direct ~ urban_dominate2 + indwork_dominate2 + v2x_polyarchy + e_migdppcln + 
                      
                      e_mipopulaln +  e_miurbani , data = data.star)
  fl.boots2[b] <- coef(boot.direct)["indwork_dominate2"]
}
sd(fl.boots)


###########################################
#######       COEFFICIENT PLOT     ########
###########################################

est <- rbind(coef(direct)["urban_dominate2"],coef(direct)["indwork_dominate2"], coef(first)["urban_dominate2"], coef(first)["indwork_dominate2"] )

#conventional standard errors
ses <- summary(first)[["coefficients"]][,2]

se <- rbind(sd(fl.boots), sd(fl.boots2), ses[2], ses[3])


#plot
pdf("CDEfig2.pdf", height=6, width=7)
plot(1, type="n", xlab="", ylab="", xlim=c(-.08, .08), ylim=c(0,10), axes=F, main="Effects conditional on foreign support to the regime")

legend("topleft", pch=c(19), cex=c(1), legend=c("CDE: IW dominate", "CDE: UM dominate", "Conventional coefficient: UM dominate", "Conventional coeffient: IW dominate"), 
       col=c("red", "green", "blue", "purple")
       
)
abline(v = 0, lty = 3, col = "black") # add verticle line

#CDE's
lines(c(est[1,] + 1*se[1,], est[1,] - 1*se[1,]), c(4,4), cex=1.2, col="orange")  
lines(c(est[2,] + 1*se[2,], est[2,] - 1*se[2,]), c(5,5), cex=1.2, col="orange")  
lines(c(est[3,] + 1*se[3,], est[3,] - 1*se[3,]), c(6,6), cex=1.2, col="orange")  
lines(c(est[4,] + 1*se[4,], est[4,] - 1*se[4,]), c(7,7), cex=1.2, col="orange")  

points(est, c(4,5,6,7), col=c("red", "green", "blue", "purple"), cex=1.3, pch=19)

# add axes and labels                                                                                         # add bottom axis
axis(side= 1)
mtext(side = 1, "Coefficient", line = 3)      

dev.off()




#####CONDITIONING ON DEFECTION

first <- lm(v2x_polyarchy_t1 ~ urban_dominate2 + indwork_dominate2 + defect + v2x_polyarchy + e_migdppcln + 
              
              e_mipopulaln +  e_miurbani ,
            data = data)
summary(first)

data$direct <- data$v2x_polyarchy_t1-coef(first)["defect"]*data$defect
summary(data$direct2)



direct <- lm(direct ~ urban_dominate2 + indwork_dominate2 + v2x_polyarchy  + e_migdppcln + 
               
               e_mipopulaln +  e_miurbani , data = data)
summary(direct)


#bootstrap the standard errors
boots <- 1000
fl.boots <- rep(NA, times = boots)
for (b in 1:boots) {
  data.star <- data[sample(1:nrow(data), replace = TRUE),]
  boot.first <- lm(v2x_polyarchy_t1 ~ urban_dominate2 + indwork_dominate2 + defect + v2x_polyarchy + e_migdppcln + 
                     
                     e_mipopulaln +  e_miurbani
                   
                   
                   
                   , data = data.star)
  boot.direct <- lm(direct ~ urban_dominate2 + indwork_dominate2 + v2x_polyarchy + e_migdppcln + 
                      
                      e_mipopulaln +  e_miurbani , data = data.star)
  fl.boots[b] <- coef(boot.direct)["urban_dominate2"]
}
sd(fl.boots)


boots <- 1000
fl.boots2 <- rep(NA, times = boots)
for (b in 1:boots) {
  data.star <- data[sample(1:nrow(data), replace = TRUE),]
  boot.first <- lm(v2x_polyarchy_t1 ~ urban_dominate2 + indwork_dominate2 + defect + v2x_polyarchy + e_migdppcln + 
                     
                     e_mipopulaln +  e_miurbani
                   
                   
                   
                   , data = data.star)
  boot.direct <- lm(direct ~ urban_dominate2 + indwork_dominate2 + v2x_polyarchy + e_migdppcln + 
                      
                      e_mipopulaln +  e_miurbani , data = data.star)
  fl.boots2[b] <- coef(boot.direct)["indwork_dominate2"]
}
sd(fl.boots)


###########################################
#######       COEFFICIENT PLOT     ########
###########################################

est <- rbind(coef(direct)["urban_dominate2"],coef(direct)["indwork_dominate2"], coef(first)["urban_dominate2"], coef(first)["indwork_dominate2"] )

#conventional standard errors
ses <- summary(first)[["coefficients"]][,2]

se <- rbind(sd(fl.boots), sd(fl.boots2), ses[2], ses[3])


#plot
pdf("CDEfig3.pdf", height=6, width=7)
plot(1, type="n", xlab="", ylab="", xlim=c(-.08, .08), ylim=c(0,10), axes=F, main="Effects conditional on military defection")

legend("topleft", pch=c(19), cex=c(1), legend=c("CDE: IW dominate", "CDE: UM dominate", "Conventional coefficient: UM dominate", "Conventional coeffient: IW dominate"), 
       col=c("red", "green", "blue", "purple")
       
)
abline(v = 0, lty = 3, col = "black") # add verticle line

#CDE's
lines(c(est[1,] + 1*se[1,], est[1,] - 1*se[1,]), c(4,4), cex=1.2, col="orange")  
lines(c(est[2,] + 1*se[2,], est[2,] - 1*se[2,]), c(5,5), cex=1.2, col="orange")  
lines(c(est[3,] + 1*se[3,], est[3,] - 1*se[3,]), c(6,6), cex=1.2, col="orange")  
lines(c(est[4,] + 1*se[4,], est[4,] - 1*se[4,]), c(7,7), cex=1.2, col="orange")  

points(est, c(4,5,6,7), col=c("red", "green", "blue", "purple"), cex=1.3, pch=19)

# add axes and labels                                                                                         # add bottom axis
axis(side= 1)
mtext(side = 1, "Coefficient", line = 3)      

dev.off()


















